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We analyze the multi-time correlations of a laser-induced 
Bose-Einstein condensate. We use quantum stochastic meth- 
ods to obtain under certain circumstances a Fokker-Planck 
equation which describes the phase-difussion process, and ob- 
tain an analytical expression of the two-time correlations. We 
perform also quantum Monte Carlo numerical simulations of 
the correlations, which are in good agreement with the pre- 
dicted analytical results. 



Therefore, multiple reabsorptions can increase the sys- 
tem energy by several recoil energies per atom [p"5[-p^. 
It is easy to understand that assuming that the reab- 
sorption cross section for trapped atoms is the same as 
in free space, i.e. ~ l/^i) the significance of reabsorp- 



I. INTRODUCTION 

During the last years the combination of laser cooling 
101 and evaporative cooling methods has lead to the 
achievement of one of the most pursued goals since the 
early days of quantum physics, i.e. the so-called Bose- 
Einstein Condensation (BEC). Such condensation is a 
direct consequence of the Bose-Einstein statistics, and 
consists in the fact that the ground state of the system 
becomes macroscopically populated. Although BEC has 
been obtained using the combination of laser and coUi- 
sional techniques, there are several experimental groups 
|-|,| which currently investigate the possibility of ob- 
taining BEC using all-optical means only. In such a case 
the number of atoms in the trap -would not decrease dur- 
ing the cooling process, and a non-destructive detection 
of BEC could be performed by simply measuring the flu- 
orescence photons. In addition, laser induced condensa- 
tion -would be easier to control externally, and could lead 
to richer effects than coUisional processes, as those em- 
ployed in evaporative and sympathetic cooling. Finally, 
the methods employed to achieve laser induced condensa- 
tion can be used to design the techniques of pumping the 
atoms into the condensate using spontaneous emission 

The most sophisticated laser cooling techniques (as Ve- 
locity Selective Coherent Population Trapping or Ra- 
man cooling are capable to cool atomic samples be- 
lo-w the photon-recoil energy, Ef( — huJn — h'^kf^/2m, 
-where /cl is the laser -wavevector, and m is the atomic 
mass. Such methods should in principle lead to BEC. 
Ho-wever, in using laser cooling to obtain the BEC, the 
reabsorption of spontaneously emitted photons turns to 
be a very important problem. This is because the sub- 
recoil laser cooling techniques are based on dark-state 
mechanisms jl4|, and the dark-states are unfortunately 
not dark respect to the spontaneously emitted photons. 



tions increases -with the dimensionality in such a 

-way that the reabsorptions should not cause any prob- 
lem in one dimension. For the case of three-dimensional 
traps several remedies have also been proposed 19 1. In 



the folio-wing we are going to assume that our system is 
one-dimensional. This assumption in the first place has 
been done for simplicity reasons, but at the same time 
it allo-ws us to avoid the reabsorbtions problem. As -we 
discuss in the conclusion section, ho-wever, our approach 
can be used in three-dimensional asymmetric traps, for 
-which heating due to reabsorbtions can also be negligible. 

When laser-cooling methods are employed, the ther- 
malization of the system occurs -with the interaction -with 
the laser field and the vacuum modes of the electro- 
magnetic field, and therefore atom-atom collisions are 
in principle not necessary. In the foUo-wing we shall con- 
sider an ideal gas (see the discussion in Sec. ||). In fact, 
the proper notion of temperature must then be revised. 
On the other hand, the dynamics of the system becomes 
more complex, because we have to deal -with an open 
system -which interacts -with the laser. In Refs. pO| , pl| , 
the quantum dynamics of a laser-cooled ideal gas has 
been studied using second quantization formalism in or- 
der to take into account the quantum-statistical charac- 
ter of the bosons, and also employing quantum stochas- 
tical methods p^|2^ to take into account the coupling 
-with the reservoir provided by the vacuum modes of the 
electromagnetic field. The equations -which describe the 
dynamics of the system have been obtained, and so has 
the stationary solution, -which corresponds for the case of 
a bosonic system -with the -well-kno-wn Bose-Einstein dis- 
tribution (BED), and under proper conditions the BEC 
appears. 

In order to characterize this laser-induced BEC, the 
f-wo-time correlations of the condensate amplitude turn 
to be very important quantities, because they describe 
the diffusion of the phase of the condensate. These cor- 
relations are very difficult to calculate both analytically 
and numerically, and until no-w no calculation of such 
correlations has been published, as far as we kno-w. The 
main aim of the present paper is to present a method to 
calculate such correlations using quasi-probability repre- 
sentations, and in particular P~representation |Q. We 
shall also present a numerical calculation based on the 
so-called -wave-function Monte Carlo method pa , -whose 
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results are in good agreement with the predicted analyt- 
ical ones. 

At this point it is worth noticing that the model con- 
sidered in this paper is not very realistic since it is one 
dimensional, and even in the relevant single dimension it 
assumes the, so called, Lamb-Dickc limit, when the size 
of the trap a is smaller than the wavelength of the laser 
field, A/27r. At the same time the model neglects not only 
elastic collisions, but also all non-elastic losses: two- and 
three-body non-elastic collisions, as well as photoassoci- 
ation losses. Realization of such situation requires modi- 
fication of atomic scattering length (as discussed in Sec- 
tion II and conclusions)and other precautions, which are 
possible (as discussed in Section II and in more detail in 
conclusions), but not easy to achieve experimentally. In 
the worst case, however, the model is realistic for small 
samples of = 15 — 50 atoms. Cooling of such small 
sample of atoms to a ground state of the trap is of fun- 
damental interest itself, and might be even useful for ap- 
plications for quantum information processing We 
stress, however, that the main achievement of this paper 
concerns methodological aspects. Even though the model 
is somewhat unrealistic, it is a paradigm model describ- 
ing quantum dynamics of the bosonic gas approaching 
the equilibrium. As we mentioned, time correlations are 
very difficult to calculate analytically and numerically in 
any model of that sort, regardless how realistic it is. We 
develop here the methods to calculate such correlations, 
and these methods are of quite importance themselves, 
since they are quite universal and can be carried over to 
more realistic models, such as those discussed by us in a 
series of Refs. [ p7| . 

The structure of the paper is as follows. In Sec. ||, 
we briefiy review the calculations already presented in 
Ref. [pOl. T he P-representation formalism is introduced 
in Sec. |II1|. The expressions obtained in Sec. |l|, allows 
us in Sec. IV to calculate the two-time correlations of 
the amplitude of the condensate. In Sec. we present 
the employed numerical method, and the comparison be- 
tween numerical and analytical results. We conclude in 
Sec. VI stressing once more methodological achievements 
of this paper, and discussing experimental feasibility of 
our model. 



II. MODEL 

In this section we briefly review the formalism already 
developed in Ref. We consider a system composed 
of N two-level identical bosons in an harmonic trap cen- 
tered at i? = 0, which are coherently driven by a stand- 
ing wave laser field, and also interact with the vacuum 
modes of the electromagnetic field. In the following we 
use, for simplicity, units with h = 1, and velocity of light 
c = 1. As in Ref. pO| , we consider four basic approxi- 
mations: (i) the laser is quasi-resonant with a particu- 
lar transition between two electronic states \g) (ground 



state) and |e) (excited state), and therefore we consider 
atoms as two-level atoms; (ii) Rotating-wave approxima- 
tion (RWA) can be used, since the frequency of the con- 
sidered electronic transition (luq) and the laser frequency 
ujl are considered much larger than any other frequency 
scale of the system; (iii) we treat the atom-field inter- 
actions in the dipole approximation since the resonant 
wavelength Ai, 3> ao, where ao is the typical atomic size; 
(iv) we neglect the atom-atom collisions, i.e. we work 
in the ideal-gas approximation. In principle, the latter 
approximation implies that we should deal with a small 
number of particles {N < 100), because we are going 
to consider that the trap is in the so-called Lamb-Dicke 
Limit (LDL), where the, so called, Lamb-Dicke param- 
eter T] = 27ra/AL < 1, with a being the size of the trap 
ground state. However, the s-wave scattering length asc 
(which governs the atom-atom interactions at low en- 
ergies) can be externally modified, for instance, using a 
magnetic field |^,^ , or a laser field |^ |, and in princi- 
ple can be made very close to zero l3l|,||, allowing 
the strict validity of the ideal gas approximation. 

We limit our discussion to the one-dimensional case. 
The relevant electronic states are then determined by the 
laser polarization perpendicular to the laser wave vector, 
which in turn can be aligned in the z-direction, — 
(0,0, wl). We assume that the atoms can only move in 
the z-direction, and are localized in transverse direction, 
so that practically one can assume that they are located 
dX R = (0,0, z). The atoms are assumed to occupy the 
energy eigenstates of the harmonic trap potential. We 
denote this energy levels by \l,g) (|m, e)) for the atoms 
in the ground (excited) electronic state, occupying the 
level l{m)— 0, 1, . . . of the harmonic trap. 

We introduce also the operators gi, gj (em, eln) that 
annihilate or create atoms in the Ith (mth) energy level 
of the ground-state (excited-state) potential. These op- 
erators fulfill standard commutation relations for bosonic 
atoms: 



(1) 
(2) 



[sm, e ,] — S„ 



Using standard quantum stochastic methods |2^-|2J] 
one can eliminate the vacuum modes of the electromag- 
netic field, and after performing a systematic expansion 
in the Lamb-Dicke parameter rj one obtains pO{ the mas- 
ter equation (ME), in the frame rotating with the laser, 
valid up to the order 0{rf'): 



p = -i[Ha + Hlas, p] + Cp, 



(3) 



where the atomic part of the Hamiltonian takes the form 

oo oo 

= ^ uj^gLgrn + ^ (wm " S)el^em, (4) 

m— m— 

with S = ujl — ujo being the laser detuning. In Eq. (^, the 
center-of-mass potentials for ground and excited atoms 
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can be well approximated by harmonic potentials of fre- 
quency = We = a; plus a small anharmonicity so that 
the energy levels are ujui + ujam'^ with to = 0, 1, . . . and 
a <C 1. From now on we use the same indices for the 
energy levels in the ground- and excited potentials. 

The Hamiltonian of the atom-laser interaction takes 
the form: 



Hlas ^V-^Yl \^"^+ 1 + eLSm+l + H.C 

rn 



(5) 



where Q, denotes the Rabi frequency describing the laser- 
atom interaction. 

The spontaneous emission part in Eq. (H) takes the 
form: 



m.m' 



_ t t 

9m' ^m'9m' 9m^m.P 



(6) 



where T is the natural line-width of the considered tran- 
sition. 

If we locate the trap at the node of the laser standing 
wave, the excitation provided by the laser then become 
weak, provided that Vl is moderate. We expect then that 
at a given instant no more than one atom will be ex- 
cited. The system is then characterized by two distinct 
time scales: the fast one, determined by NT, uj and 5\ and 
the slow one, characterized by TfVt'^/NT,rfVt^/{uj ± 5), 
which describes the jumps between the various ground- 
state levels. In such a case one can use standard adiabatic 
elimination techniques ||2^ to remove excited state popu- 
lations. Physically, this reflects the fact that a combined 
process of excitation and (relatively rapid) spontaneous 
decay causes a redistribution of atomic population among 
different levels of the ground-state trap. The ME after 
the adiabatic elimination of the excited states becomes: 



(7) 



where the coupling between the levels m and to -I- 1 ||3^ ] 
is given by [g^ : 

+ ^("1 + 1)(2ALM™ - AraAl^p - pA„^Al), (8) 



with Ajn = gla9'm+i, and 



(iVr/2)2 + (a. T Sf 



(9) 



It is easy to check that for (5 < there exists 
an exact steady state solution of ME (|^) which has the 
canonical form: 



Pst = -^9^'" 



(10) 



where q = exp(— w/fcsT) = r+/r— , and Z is the canon- 
ical partition function. For the particular case of a one- 
dimensional harmonic potential the partition function 
can be derived in closed form fSSll : 



^ 1 



(11) 



The thermodynamics of the system is determined by the 
Helmholtz free energy F = —kBT\n{Z). The canonical 
version of the chemical potential is then given by: 



M = idF/dN)T,^ = fcBTln[l - 
There is a non-zero temperature [^6|| 
hcj N 



kB HN) ' 



(12) 



(13) 



below which p is effectively zero, q can be conveniently 
re-written in terms of this critical temperature: 



q = exp 



HN)_Tc 
N T 



(14) 



III. P-REPRESENTATION 

In the following section, we are going to employ ME 
(0) to calculate the two-time correlations of the ground 
state in the stationary regime, i.e. correlations of the 
form {gQ{t)go{0)). In order to do that it is convenient to 
use the so-called Glauber-Sudarshan P-representation 
p3[ , which is defined as: 

p{t) = J (fzQ(fzi . . .P{zo, zI\zy,z\\ . . .) 

X \za,zi,...){zQ,zi,...\. (15) 

where \zq, zi, . . .) are coherent states. After transforming 
into P-representation, Eq. (0) becomes: 



P = 



OO 



(16) 



m=0 



where the transition between to and to -t- 1 is given by 
the operator: 



Sm = — (to + 1) 



d 



Zm+l I 



dZm+l 



+ — 

' dzt 

Zm+l 



\Zm\' 



+ 1 



\Zm+l\ 



dz„ 



3 



2q 



\Zm+l\ 

\Zm+l\ 
92 



m+l 

d 

dZm+l 



Zm + 1 + 



1 + 



- \Zrn\ - 



d 
dz„ 



dz* 



Zm 1 |Zm+l| — 



Zm+l\ 



d 



d 



dz, 



Zm+l + 1 
'Zju+l + 1 



m+1 



(17) 



Let us define the reduced P-representation for tlie 
ground state of the trap as: 

Po(zo,Zo;0 = / d?ziSz2...P{zQ,zQ*]zl,zl*]...]t). 



(18) 

The equation which describes the dynamics of Pq, takes 
the form: 

Po = —Q { —^0 + Q^xo ) Po 



+ ^(9-1) 



d 



2 ' \dxa 

N dxodxl 



xo 



^0 Pi 



(19) 



where Pi — j (Pzi<fiz2 ■ ■ ■\zi\'^P. We have also intro- 
duced here the convenient notation zo ~ \/Nxo. 

One can prove (see Appendix ^) that for times t > 
[r_(l — q)N\xo\'^]~-^ , one can adiabaticaUy ehminate the 
excited trap states, and retrieve a closed Fokker-Planck 
equation (FPE) for the reduced Pq: 



Po = 
1 

~4 



N{1 — q) { dxodxQ 



d 



d 



"5 — ^0 + ^^^0 

OXq OXq 



dx( 



-Xq 



dXn 



Fol 



Po- (20) 



Changing into the more convenient variables zq 
rexp{i(jji), the FPE (|2C|) becomes: 



P0--—P0, 

j,2 QM2 



with 



4(1 -g)- 



(21) 



(22) 



Eq. ( |2l| ) is one of the central results of this paper, 
and has a simple physical interpretation as an equation 
describing diffusion of the phase of the condensate wave 
function. Interestingly, even though we are dealing here 
with an ideal gas approaching the thermal equilibrium, 



the equilibrium state has temporal properties character- 
istic for the states obtained via spontaneous breaking of 
the U{\) phase symmetry, such as in the theory of laser, 
or in general any theory of second order phase transi- 
tions in which the effective potential for the order pa- 
rameter below the transition temperature has a Mexican 
hat shape |Q. In particular, the phase diffusion rate 
is here inverse proportional to r^, i.e. to the number of 
atoms in the condensate. 



IV. TWO-TIME CORRELATIONS OF THE 
CONDENSATE 

The solution of the FPE ( ^l]) can be obtained using 
the following Green function: 



G(r,0,i||r„,(/.o,i = O) 



E—n -^t 
cosn(0 — 0o)e 



(23) 



n=l 



Using this we can already calculate the time correlation: 

{9l\t)9om 



27r 

rdrd(p 
Jo 



2tt 



rodrod(l)o 

/o "'0 "'0 ^0 

xG(r, i||r„, 00, t - 0)Po(ro, 0o)r-V^e-^'('^-'^'') 



rdrdcjyr'^'' e ' ^*Po(r, 



(24) 



where ":" denotes normal ordering. In the following we 
just consider the case of an exponent 1 = 1. We expand 
the exponent, assuming that for temperatures T suffi- 
ciently below Tc, fluctuations of no = 9o9o sre small: 



la 

no 



ICt 



(no) - ((n-o) - no) 
ICt ( ^ ( ^ no 



(no) 
2Kt 



1 



Kt 



(no) {noY 



1 - 



(no) 



(25) 



Using this expression and the identity : exp(^(7j(7o) ■ = 
exp(ln(^-)- l)(7j(7o)i Eq. ([2^) can be rewritten in the form: 



(5i(i)ffo(0)) 
= (: 5^50 exp 



2/Ct 



9I90 



[no) [no) 

^-2Kt/(no) _ 1) 1^(1 + Ja/{noY)\ ) . (26) 



Note that from Eq. (^6|), it is clear that the correlations 
depend on the fluctuations of the condensate fraction, 
which are not correctly described using the usual text- 
book treatment of the system based on a grand-canonical 



4 



ensemble 35-3^. In order to calculate the correlations, 
it is thus necessary to calculate the required averages 
using physically sound ensemble, which in this case is 
the canonical one. Using canonical ensemble we have to 
determine the probabilities to have hq particles in the 
ground state, given a total number N of particles [EgI: 



N 

n (1 

j=N-no + l 



(27) 



Therefore: 

{(glmgom) 



N 



^-2Kt/{no) 



Po^^(no|A^)no 1 



"0=0 



_ja_ 

{no)' 



rio-l 



(28) 



The above closed analytic formula is another central re- 
sult of this paper. 



V. NUMERICAL RESULTS 

Numerical simulations of the system under considera- 
tion can be performed using Quantum Monte Carlo tech- 
niques. In our case, it is a quite difficult task due to the 
large size of the Hilbert space which has to be simulated, 
even if only a relatively small total particle number N 
and a low cutoff of the trap level structure is chosen. 

To avoid having to compute the dynamics of the huge 
density matrix of this system, we applied a stochastic 
wave function method, which replaces the deterministic 
evolution of the density matrix (following Eq.(|^)) by a 
stochastic evolution of an ensemble of state vectors pO| . 

A naive application of this algorithm to simulate (|7|) 
would lead to a stochastic jump process in which the op- 
erators Ajn and A]^ act as jump operators. This process 
would represent the density matrix as an ensemble of si- 
multaneous eigenvectors of all occupation number opera- 
tors n,n — glndrm Causing large fluctuations of the phase 
which we are interested in. Consequently, the Monte 
Carlo simulations would require a very large number of 
realizations. 

We found that convergence can be improved by rewrit- 
ing the Liouvillian operators given in Eq. (^) in the fol- 
lowing way: 



with operators Am'^ := -^[Am ± ?7m), where the r]m are 
c-numbers which are chosen to be of the same order of 
magnitude as the mean occupation number -y/ {rimnm+i) , 
so that rjm ~ Am- Eq. ( p9| ) has an advantage that while 
it is, of course, reproducing the same master equation 
Eq. (^, it leads to a different stochastic jump process, 
in which the Am^ act as jump operators. This process 
represents p as an ensemble of vectors which are, in gen- 
eral, superpositions of different occupation number eigen- 
states, and fluctuations of g{t) within this ensemble are 
much smaller. This observation is yet another important 
result of this paper. 

Due to computational constraints we have simulated 
the case of iV — 10 particles, confined to the 5 lowest 
levels of the harmonic trap, for different values of the 
temperature parameter q = T^/T^. For each value of q, 
the function g{t) was estimated by averaging over 1500 
trajectories. 

Figure (|l|) shows the results of the numerical simula- 
tion, compared to the analytic formula (p8|), for differ- 
ent values of q. Quite remarkably, the analytical result 
(based in the approximation of (no) ^ N — (no)) is in 
very good agreement with the numerical results even for 
N equal just 10 particles (for larger number N the agree- 
ment should be even better). As expected, the agreement 
is much better for low temperatures than for high temper- 
atures, where the numerical simulation indicates a faster 
decay of g{t) than that predicted by the low-temperature 
approximation (|2F 



2 ^ 


n){2Ai+^pAW^ 




-paw'aW) 


r_ , 

+— (m4 


-l){2A(JpAi-^^ - 








-1)(2AWV4,I' - 


rn m r 


pA^^A^i^') 




-1)(24(-)V^L^ - 


m m r 


pA^jA^m^^), (29) 



VI. CONCLUSIONS 

In this paper we have studied quantum dynamics of 
a Bose gas in a trap undergoing sideband laser cooling 
in the Lamb-Dicke limit. The master equation describ- 
ing the dynamics of the system can be regarded as a 
paradigm equation for collective cooling dynamics. One 
of the difficult problems associated with such dynamics 
consists in calculating time dependent correlation func- 
tions, such as for instance those that describe temporal 
phase fluctuations of the Bose condensate. In this paper 
we have presented a solution to this problem. Our main 
results should be regarded from the methodological point 
of view: 

• We have formulated an analytic method to describe 
temporal correlations based on an expansion of the 
master equation valid at low temperatures, when a 
large number iVp of particles are in the condensate. 
The expansion parameter in our approach is I/Nq. 
The expansion can be, and already has been applied 
for other models of the Bose gas dynamics, that 
describe more realistic physical situations (see Refs. 

We have formulated a numerical method to calcu- 
late the time correlations which modifies the jump 
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processes involved in the master equation in such a 
way that the corresponding Quantum Monte Carlo 
simulations are much more stable and require av- 
eraging over much less number of quantum trajec- 
tories to achieve good accuracy. The method pro- 
posed is general, and can be applied not only for 
the present model, but also for more realistic re- 
lated models. 

• Analytic and numerical results agree very well even 
at the border of the validity of the analytic theory 
{No ~ 10). 

• Collective laser cooling leads at low temperature 
to phase diffusion in a Mexican hat potential. This 
result is also general, and holds for other more re- 
alistic models of the collective laser cooling. 

The above listed methodological results provide the 
main value of the present paper. It is, nevertheless, inter- 
esting to speculate whether the consider model is purely 
academic, or whether it can be realized experimentally. 
We shall argue now that it can be realized for small num- 
ber (iV = 10 - 50)) of particles. 

Let us discuss step by step the most relevant approxi- 
mations used in this paper. 

• 1 Dimension. This approximation was done mainly 
for technical reasons. The model can be easily 
generalized to describe condensation in a three- 
dimensional trap. In that case, additional precau- 
tions should be taken into account to avoid the re- 
absorbtion problem. In particular, if the width of 
the excited states is smaller than fiui, we were thus 
working not only in the LDL limit, but also in the 
festina lente regime ]l9| ]. Extending the dynamics 
to three dimensions in an asymmetric trap will in 
this regime not introduce any kind of dangerous 
reabsorbtion problems. In fact, the dynamics will 
consist of 3 independent dynamics corresponding 
to the cooling in the x, y, and z direction. 

• Absence of elastic collisions. As we have men- 
tioned, this requires that the atomic density should 
be sufficiently small, or alternatively that the scat- 
tering length should be modified to very low val- 
ues. The first possibility is not interesting, because 
we require also that the LDL conditions are ful- 
filled, which means that a should be of the order 
of at most O.l^m. The condition to be fulfilled is 
AirNh'^asc/ fnV < with the effective conden- 
sate volume V — (27r)^/^a'^. For the parameter 
considered this gives gives N( < 20, where ( is the 
modification factor of the scattering length. 

• Absence of two body inelastic collisions This prob- 
lem has a simple remedy. The cooling and con- 
densation should take place in a dipole trap, and 
occur in the electronic ground state of the atoms. 



Two body inelastic processes are then completely 
suppressed. 

• Absence of three body inelastic processes Three 
body losses can be typically neglected provided the 
density is less than 10^^/cm^. For a = 0.1/im that 
requires N < (27r)^/^ ~ 15. If the three body losses 
are modified in a corresponding way as elastic col- 
lisions (which seems to be the case for *'^Rb p2|), 
then the corresponding condition is much less re- 
strictive NC^ < 15. 

• Absence of photoassociation losses This can be re- 
duced by using red detuned laser tuned in between 
the molecular resonances ||3^ . Even is such a case, 
photoassociation losses become relevant when the 
density reaches the limit lO^^/cm"^, i.e. in our case 
for N < (27r)^/^ ~ 15. It seems likely, that this es- 
timate can be improved significantly when the laser 
is tuned below the minimum of the molecular tran- 
sition. While it can hardly be though of for the 
direct transition, we should stress that laser tran- 
sition considered in this paper can be equally well 
regarded as a Raman transition. In such case, tun- 
ing of the stimulated two-photon transition below 
the minimum of the molecular resonance is possi- 
ble. 

Summarizing, we see that even in the worst case (no 
modifications of the scattering length, no special precau- 
tions regarding photoassociation losses) our model should 
be valid for ~ 10 — 15 particles. Additional precau- 
tions can allow to extent the validity of the model to 
significantly larger values of N. 
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APPENDIX A: ADIABATIC ELIMINATION OF 
THE NON-CONDENSED STATES 

In this Appendix we present in detail the calculations 
which allow us to transform Eq. ( p^ ) into Eq. (pO|). First, 
we show that our problem can be reduced to a two-level 
system formed by the level and 1 of the trap. Then, we 
adiabatically eliminate the level 1. 



2r_ [(1 - q){ni)Gc{n2)GC - qMgc + Mac] = 0. 

(A7) 

Therefore the contribution of iSi for Pi cancels out. In 
general, the contribution of Si for P" is not exactly zero, 
but it is always a constant, following the MD arguments. 
Therefore our system reduces to a two-level system, in 
which only the levels and 1 must be considered. 



1. Two— level system 



2. Adiabatic elimination of the level 1 



First, we shall analyze the dynamics of Pi . From jl^ ) 
and (0) one obtains that 



Pi = 



N 



V^ + (g-l)/o 



2(<7-l)iV|xo|2 + (2g+l)?o-2 



Pi 



-T^qN\xo\^Po 



-2r_ 



(1 - 9)^1,2 - iPl + P2 



where we have used the notation: 
92 



dxodxQ ' 
f d d , 



(Al) 

(A2) 
(A3) 
(A4) 



and Pi^ 



p2,0 



pO,l 
^1,2 ■ 



Let us analyze in detail 



the last line in the RHS of Eq. (Al), which comes from 
the contributions of Si in Eq. (16), i.e. the contribu- 
tions given by the transitions 1 <-> 2. For a tempera- 
ture sufficiently below Tc, and sufhciently large N, the 
atoms in the non-condensed states of the trap form a 
so-called Maxwell-Demon (MD) ensemble |3^], i.e. an 
ensemble which exchanges particles with a reservoir pro- 
vided by the condensate without exchanging the energy. 
Therefore, the excited states can be considered as: (i) 
independent of the population of the ground state; (ii) 
decorrelated among each other. Due to these properties: 



Pi 



.j.k 



1,2 



{ni)Gc{n2)GC, 



(A5) 



where the subindex GC means that the averages are cal- 
culated in the grand canonical ensemble. It is well known 
that these averages have the simple form: 



Having reduced the system into just two levels, and 
1, we shall adiabatically eliminate the level 1. Thi s can 
be achieved because, as observed from Eqs. ( |AlD and 
(|l9|). Pi decays on a time scale of the order 0{1/N), 
whereas Pq decays in a time scale of the order 0(1). We 
are interested in contributions up to the order 0{1/N) in 
Pq, and hence in contributions up to the order 0{1/N) 
in the stationary value of Pi. From Eq. (|Al| ) it is clear 
that the contributions of the order 0(1) in Pi lead to 
terms of the order 0{1/N) in the stationary value of Pi. 
Therefore we are interested in Pi up to order 0(1), and 
therefore in Pi^ up to order 0(1). This means that we 
need to calculate Pf just up to the order 0{N): 



P? = [-Ml - q)N\xon Pi 
+ I_ [8qN\xo\^] Pi. 



(A8) 



For times t > [2r_(l - q)N\xQ\'^]~^ , the stationary val- 
ues: 



Pi = 2 



Pi- 



rn 



is obtained. Therefore Eq. ( [Al[ ) reduces to the form (up 
to order 0(1)): 



-2{1 - q)N\xo\^ + lo - 2 



Pi 



+ — [2qN\xo\^] Po. 
And for t > [r_(l - q)N\xo\^]-'^, 



(AlO) 



Pi 



Po- 



l-q 2A^(l-g)2 °|xoP 

Substituting this value in Eq. (p^), one obtains Eq. (|2 



(All) 



j/GC 



l-qi 



(A6) 



Therefore, the last line in the RHS side of Eq. (Al) be- 
comes: 



8 



50.0 10D.0 150.0 200.0 



(b) 



25.0 50.0 75.0 100.0 



(c) 



25.D 50.0 75.0 100.0 



FIG. 1. Comparison between the analytic approximation 
described by Eq. ( [28[ ) (dashed line), and the results of 
the numerical simulations (solid line), for different values of 
q = r+ /r_ . In the simulations we have chosen the total par- 
ticle number A*' = 10. The time scale is in units of TZ^- The 
numerical results have been obtained by averaging over 1500 
trajectories for each value of q. 



9 



